Universal Thermometry for Quantum Simulation 
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Quantum simulation is a highly ambitious program in cold atom research currently being pursued 
in laboratories worldwide. The goal is to use cold atoms in optical lattice to simulate models for 
unsolved strongly correlated systems, so as to deduce their properties directly from experimental 
data. An important step in this effort is to determine the temperature of the system, which is 
essential for deducing all thermodynamic functions. This step, however, remains difficult for lattice 
systems at the moment. Here, we propose a method based on a generalized fluctuation-dissipation 
theorem, ft does not reply on numerical simulations and is a universal thermometry for all quantum 
gases systems including mixtures and spinor gases. It is also unaffected by photon shot noise. 



At present, there is worldwide experimental effort to 
simulate theoretical models for strongly correlated quan- 
tum systems using cold atoms in optical lattices. If suc- 
cessful, these simulations can provide detailed thermo- 
dynamic information for many models whose solutions 
are unknown, even some of them (such as 2D Hubbard 
model) have been studied for decades. To deduce the 
thermodynamic properties of these models directly from 
experiments, it is necessary to determine three quantities 
accurately : density n, chemical potential /i, and tem- 
perature t m m. The recent experiment of Cheng Chin's 
group Q using in situ density profile to identify directly 
the thermodynamic phases for boson Hubbard systems is 
a very important step toward realizing the full power of 
quantum simulation [lj. The prospect of this realization 
is further enhanced by the impressive improvement in res- 
olution of density imaging recently developed in Markus 
Griencr's group[4j. The next crucial step is to have an 
accurate temperature determination. 

Often, the temperature of a lattice gas is estimated 
by assuming the lattice is turned on adiabatically. One 
then equate the entropy of the final state S/(Tf) to that 
of the initial state Si(Tj), (i.e. the state before the lattice 
is switched on), and then deduce the final temperature 
Tf from the initial temperature Tj through this relation. 
One factor detrimental to this procedure is the intrinsic 
heating caused by spontaneous emission, which occurs as 
the lattice is turned on, and during the time when exper- 
iment is performed Q. To make things worst, the entropy 
function Sf(T) of many systems of interest remains un- 
known. So the errors of this method are uncontrolled @. 

For quantum gases in a single trap without optical lat- 
tice, their temperatures can be deduced from the density 
profile at the surface, which has the Boltzmann form. 
In principle, one can apply the same method for lattice 
quantum gases, as interaction effects becomes unimpor- 
tant near the surface. However, an accurate determina- 
tion of the density profile near the surface will require 
improving the imaging resolution to a single site. It 
will also require repeating the experiment many times 
so as to achieve a good signal to noise ratio. To avoid 
these demands, many experiments resort to the afore- 



mentioned adiabatic assumption for temperature deter- 
mination. However, due to the uncontrolled errors in 
this method, it is desirable to have an alternative scheme 
which is robust and free of all the problems mentioned 
above. We also note that by studying the density at the 
surface, one can not determine whether the entire sample 
is in global equilibrium. 

In this paper, we present a new scheme to determine 
the temperature of trapped quantum gases based on the 
fluctuation-dissipation theorem for non-uniform systems. 
This method applies to all quantum gas systems (single 
component gases, mixtures, spinor gases) and is unaf- 
fected by background photon shot noise. It can also be 
used to deduce magnetic susceptibility of bulk systems. 
This method does not require numerical input, and can 
tell whether the system is in global equilibrium. 

Al. The proposal: We begin with two basic assump- 
tions used in most experiments on quantum gases which 
have been justified in many cases. The first is that the 
density n(r) of a quantum gas in a trap V'(r) can be calcu- 
lated in grand canonical ensemble , i.e. n(r) = n(r; T, fi), 
where 



n(r;7» = 



Trri(r)e" 



Tie -0(H+V-nN) 



<n(r)> 



7>- 



(1) 



where (3 = l/(ksT), H is the Hamiltonian without trap- 
ping potential, T is the temperature and fi is the chemical 
potential. The second is that n(r; T, //) is given accu- 
rately by local density approximation (LDA), i.e. 

n(r;T,/i)=n (/i(r),T) ) M (r) = /, - V(r), (2) 

where n D [y, T) is the density of a homogeneous sys- 
tem with hamiltonian H and chemical potential v, i.e. 
n (y,T) = Tre-^-^N/iQTre-^ 6 and ft is 
the volume of the homogenous system. For lattice quan- 
tum gases, LDA is justified if the variation of trapping 
potential between neighboring sites is small compared 
with the hopping matrix element. Eq.JT]) implies 



k B T 



d(n(r)) 
<9/i 



dr'[(n(r)n(r'))-(n(r))(n(r'))], (3) 
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where 
V(v) = 



<->t>. 
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For an isotropic harmonic trap 



^Muj 2 r 2 with frequency u, Eq.© becomes 



k B T d(h(r)) 
Mia 2 r dr 



= /dr'[(n(rHr'))-(n(r))(n(r'))], 

(4) 



or simply 



k B T 0(n(r)) 
Mu! 2 r dr 



(n(r)N) - (n(r))(N). (5) 



Eq. ([5]) suggests a convenient way to determine temper- 
ature. Suppose we repeat the experiments Q times, and 
label the measured quantities of each sample by a su- 
perscript "i", i = 1,2,3, ...Q. Let nW(r) be the density 
profile of the i-th sample, and — f (r) be the 
total number of particle of that sample. The averages 
of these quantities over all Q samples will be denoted as 
n(r) and N, where x = J2i=i x ^ IQ- l n * ne limit where 
Q >> 1, Eq.© can be written as L(r) = R(r), where 
R(r) = n(r)N - n(r) N, or 



R(r) =Q- 1 J2n^ l \r)N^ +Q- 



Q 

n^{r)N^; (6) 
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dr 
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That we label temperature T with a superscript i is be- 
cause in real experiments, there are fluctuations in tem- 
perature in each of the samples due to the initial evap- 
oration process. In the following, we shall assume the 
temperature fluctuations from sample to sample are suf- 
ficiently small compared to the mean temperature so that 
they can be ignored. In this case, we can set to its 
mean values, which we simply denote as T, and Eq.© 
becomes 



L(r) = T£(r), £(r) 



k B dn(r) 
Muj 2 r dr 



(8) 



Eq.© them implies T = R(r)/£(r) at any position r. 

There is, of course, the practical matter of how many 
samples is needed to average over to reach the thermal 
average. A very large value of Q will not be practical. To 
achieve fast convergence, one can suppress the noise by 
averaging over a ring of thickness e, resulting in the func- 
tion (say, in the 2D case) ((p) — f p+t r](p')p'dp' /Q(p), 

where rj(p) — d9 n(p, 9) is the angle integrated den- 
sity at radius p, Cl(p) = n[(p + e) 2 — p 2 ] is the area of the 
ring being averaged over, and (p, 9) are polar coordinates. 
From Eq.©, we obtain temperature T as 



T = K(p)/£(p), 



K(p) = ((p)N - C(p) N, 



(9) 



(10) 
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FIG. 1: Density profile in the trap. Red dots: Ensemble aver- 
aging over 2000 configurations. Blue diamonds:LDA. Purple 
boxes: Exact density in the trap. 



Mcj 2 n(p) 
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Eq.© holds for all radii p. In the 3D case, the 
quantity easily accessible is column integrated density. 
The corresponding expression for Eq.© is to replace 
C,(p) and r](p) by their column integrated analogs a(p) 
and b(p), where a(p) = J p p+t p'b{p')dp' /n{p), b(p) = 

d9 J dz n(p,9,z), and (p, 9, z) are cylindrical coor- 
dinates. 

To illustrate the working of Eq.©, we consider a 2D 
ideal Fermi gas in a square lattice with Hamiltonian 



t 



and in an overall harmonic 



potential V = |Er^ w2 ^ 2c r c r with frequency w. 
Here, t is the hopping matrix element, (R, R') means 
neighboring sites, and cj^ creates a fermion at site 
R with spin a. The equilibrium density of this non- 
uniform system is (n(r)) = J2 a \ u a{ r )\ 2 f{E a ), where 
f{x) = {e\ x ~i J ->' kBT + 1) _1 is the Fermi distribution func- 
tion, E a and u a (r) are eigen-energies and eigen- functions 
of the system H + V. 

In Figure 1, we show the equilibrium density of a sys- 
tem with temperature T/t = 0.1 and a chemical potential 
/i adjusted so that the number of particles is N — 1200. 
We also show on the same plot the LDA result, which 
differs from the grand canonical result by less than 0.1% 
and is invisible in the figure, justifying the assumptions 
we mentioned at the beginning. To generate an equilib- 
rium ensemble, we start with an arbitrary assignment of 
and 1 of the occupation numbers {n a } of the energy 
levels {E a } up to a very large cutoff A, and evolve the 
set {n a } with Monte Carlo scheme. We have generated 
2000 configurations {n a } after the system has reached 
equilibrium, which we refer to as the equilibrium ensem- 
ble. The average of the occupation number n a is given 
by the Fermi distribution (with 0.1% accuracy in our 
calculation), and that there are no correlations between 
the occupations of different energy levels (which is the 
properties of ideal gas). Within this equilibrium ensem- 
ble, we randomly select Q configurations, which would 
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FIG. 2: Compressibility (left) and number fluctuations (right) 
in the trap. Red dots: averaging over 50 configurations. Pur- 
ple boxes: Exact results in the trap. Blue diamonds: LDA. 



correspond to Q measured samples in experiments. The 
level occupation of these configurations will be labeled 
as {ria }, i = 1 to Q. We then evaluated the densi- 
ties n(r) = J2 a \ u a(r)\ 2 na and the number fluctuations 
n(r)N — n(r) N — ^2 a |u a (r)| 2 (n 2 — n^ 2 ) of these Q sam- 
ples, and the angle averaged quantities lZ(p) and C(p) de- 
fined in Eq.lUnD and JITJ, where ^ = Q^ 1 Ei=? n ^ and 



q- x e:=? 



Q^M 2 



The dependence of these quanti- 



ties as a function of raduis p for the case of Q = 50 is 
shown in Figure 20. In this figure, we also include the 
LDA result, which differ from the data only by a few 
percent. This shows once again that LDA is an excellent 
approximation. Figure 3 displays the pairs (1Z(p) , C(p)) 
in the C — 7Z plot. According to Eq.©, all the points 
should fall in a straight line with slope given by T. By 
randomly choosing 50 samples in the equilibrium ensem- 
ble, we obtain a temperature within 3% of the actual 
temperature. If we increase the number of sample to be 
averaged to Q = 200, the accuracy in temperature in- 
creases to 1%. We have repeated our calculation for the 
same system at lower temperature T/t = 0.02 and have 
found the same accuracy in temperature determination. 

We stress that the angular average is crucial for our 
scheme. Due to the self-averaging property of the equi- 
librium ensemble, the angular average enhances signal to 
noise significantly, and amounts to a significant increase 
in the number of configurations averaged. The high ac- 
curacy of temperature determination by averaging only 
50 samples makes our scheme practical. We would also 
like to point out that if the system is not in global equi- 
librium, but was able to establish different temperatures 
in different parts of the sample, then the points in the 
1Z-C plot will fall into a few straight lines with different 
slopes. 

A2. Fluctuations in p and T: Next, we consider 
the effect of fluctuations in T. We have generated den- 
sity profiles of equilibrium ensembles at different tem- 
peratures while keeping p and lu fixed. We find that 
with 1% (5%) temperature fluctuations, the accuracy for 
temperature determination after averaging 50 samples 
remains at 5% (changed to between 5% to 10%). We 
have repeated our calculations for similar fluctuations 
in p (which amounts up to 5% fluctuations in N), and 



have found the results. This shows our scheme is robust 
against these fluctuations, and that Eq.© is justified. 

A3. Local density fluctuation : We again consider 
the 2D case. In experimental analysis, one divides up the 
real space into a square lattice of units cells (i.e. bins) and 
count particle number n R where R labels the location of 
the bin. Eq.([3]) then becomes 



R 



D 



R 



R- 



(12) 



where Cr and -Dr are the local compressibility and local 
density fluctuation, 



C R = - 



Muj 2 r dr 



Dn = (n R ) - <n R > 2 , (13) 



Fr = Er/-^r Gr,r' is the fluctuation due to neighboring 
bins, and Gr.r< = (n R n RI ) — ( n R )( n R ') is the density 
correlation at different bins. If the range of Gr.r' at R 
happens to be very short, which may occur if the system 
at R is in the Mott phase or the normal phase, then 
Fr ~ and Eq. (fl"2"|) implies that temperature is simply 
the ratio T = D R /C R . Thus, if we plot D R = J e D R 
against Cr = J e C R , where J g denotes angular average, 
we will find many points (Cr, Dr) in the C — D plot fall 
onto a straight line while many other points do not. The 
former comes from the regions of {R} with short range 
density correlations, while the latter from regions with 
longer range correlations. This suggests a simple way to 
use local density fluctuation to determine temperature: 
Even though only a portion of the curve (Cr, Dr) fall on 
a straight line, one can still determine T from the slope 
of this straight line[8|. This method, while convenient, is 
not as general as that discussed in A.l, for it requires a 
significant part of the system to have short range density 
correlations. 

A4. A further simplification: Finally, we note that 
by integrating overall r, and using the fact that dp = 
-±Mw 2 dr 2 , Eq.© becomes (in the 2D case) 



2irkBTn R= 
Muj 2 



(n 2 ) - (n^ 



(14) 



One can therefore also determine T from the central den- 
sity and total number fluctuation. In the 3D case, n R= o 
will be replaced by the column density at R = 0. De- 
spite the simplicity of Eq. (fT4"]) . the algorithm discussed 
in Al remains more robust, as T is determined by the 
contributions of all radii p. 

B. Photon shot noise: In real imaging process, there 
is photon shot noise. The measured atom number (de- 
noted as n ex (r)) is related to the the actual atom num- 
ber h(r) as n ex (r) — n(r) + v(r), where v(r) is the con- 
tribution due to photon shot noise. Since photon shot 
noise is a property of the laser, its probability distribu- 
tion V p [v\ is independent from that of the density distri- 
bution, V a [n], which is given by thermodynamics. Av- 
eraging over different equilibrium samples (denoted as 
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FIG. 3: Linear fit for {£(p),lZ(p)} to extract T. Blue dots are 
the results of C(p) and TZ(p) from averaging 50 configurations. 
Blue straight line is the fitting results. Red straight dashed 
line represents the real temperature T = O.lt. 



(~)a,p) means averaging over these two independent dis- 
tributions. In other words, the experimentally measured 
density n ex (r) = (h ex (r)) a ^ p is 

n ex (r) = n(r) + c, n(r) = (n(r)>„, c = (v{v)) p , (15) 

where c the average of background photon shot noise, 
which can be calibrated by taking images in the absence 
of atoms. We shall also assume the short noise has no 
significant spatial correlation, i.e. 



{v{v)v(r')) a = c 2 + ci<5(r -r'), 



(16) 



where C\ is shot noise fluctuation about its mean c. 
Since the average noise c is independent of /i, we have 

fcsT 9n-(r) = fcsT 5n(r) = {f . ^ _ ^ 



N = (N) a = J drn(r). Next, consider the fluctu- 
ations of measured density, i? ex (r) = (n ex (r)N ex ) a<p 
-{h ex (r)) a , p {N ex ) a , p . We note that 



(h ex (r)N ex ) a , p = ((n(r) + P(r))(N + J dr'i>(r'))>«,p(18) 
= (h(r)N) a +cN + n{v)Vc + J dr' ' {u(v)u{v'))) p (19) 

where V is the volume for photo collection; and 

(n ex (r)) a , p (N ex ) a , p = (n(r) + c)(N + Vc) (20) 
= n(r)N + cN + n(r)Vc + Vc 2 . (21) 

Subtracting Eq. (f21~j) from (fl~9|| . and using Eq. (fl6|) . we 

have 

R ex (r) = {n(r)N) a -n(r)N + Cl . (22) 
Eq.(|T7I) and (J22J) then imply 



dn ex (r) 
9/i 



(h ex (r)N ex ) a . p - (n ex (r)) a , p (N ex ) a , p - Cl . 



T 



\ k B ) rP+* 



C ex (p) N - ci 



J p p+e dsdr] ex (s)/ds 



(24) 



(p) and Ci shall be replaced by 
z)/Cl(ft) and 



In 3D 

a ex (p) = f p P+e dp'p>J™d6jdzn™(p>, 
b ex (p) = J w d6 J dz n ex (p, 6, z) and C\ — f dz C\ respec 
tively. 



Conclusion: We have shown that density fluctuation 
is a powerful way to determine the temperature of a 
trapped gas. It is clear from our derivation that this 
method applies to other systems such as mixtures and 
spinor gases. The fact that the temperature can be de- 
termined by the fluctuation at every point in the sam- 
ple provides considerable cross checks on the accuracy of 
the result. Our method can also reveal situations where 
different regions of the sample are in equilibrium within 
themselves but not with each other. At present, all meth- 
ods of thermometry requires the input of specific theo- 
retical modeling. Our method replies only on thermody- 
namics. It is therefore immune from errors of theoretical 
modeling, and is in line with the true spirit of quantum 
simulation, i.e. finding information of unsolved models 
without specific theoretical input. 
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